{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "combined-semiconductor",
   "metadata": {},
   "source": [
    "# The Freshman Plague"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-table",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "electoral-turkey",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# install Pint if necessary\n",
    "\n",
    "try:\n",
    "    import pint\n",
    "except ImportError:\n",
    "    !pip install pint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "formal-context",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://raw.githubusercontent.com/AllenDowney/' +\n",
    "         'ModSimPy/master/modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "progressive-typing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "undefined-miller",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/' +\n",
    "         'chap11.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "growing-sperm",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import code from previous notebooks\n",
    "\n",
    "from chap11 import make_system\n",
    "from chap11 import update_func\n",
    "from chap11 import run_simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "plastic-trigger",
   "metadata": {},
   "source": [
    "[Click here to run this case study on Colab](https://colab.research.google.com/github/AllenDowney/ModSimPy/blob/master/examples/plague.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "complex-renewal",
   "metadata": {},
   "source": [
    "This case study picks up where Chapter 12 leaves off."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "recent-cooper",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_immunization(system, fraction):\n",
    "    system.init.s -= fraction\n",
    "    system.init.r += fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "found-learning",
   "metadata": {},
   "outputs": [],
   "source": [
    "tc = 3             # time between contacts in days \n",
    "tr = 4             # recovery time in days\n",
    "\n",
    "beta = 1 / tc      # contact rate in per day\n",
    "gamma = 1 / tr     # recovery rate in per day\n",
    "\n",
    "system = make_system(beta, gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "synthetic-element",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_total_infected(results, system):\n",
    "    s_0 = results.s[0]\n",
    "    s_end = results.s[system.t_end]\n",
    "    return s_0 - s_end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "changed-capability",
   "metadata": {},
   "source": [
    "## Hand washing\n",
    "\n",
    "Suppose you are the Dean of Student Life, and you have a budget of just \\$1200 to combat the Freshman Plague. You have two options for spending this money:\n",
    "\n",
    "1.  You can pay for vaccinations, at a rate of \\$100 per dose.\n",
    "\n",
    "2.  You can spend money on a campaign to remind students to wash hands\n",
    "    frequently.\n",
    "\n",
    "We have already seen how we can model the effect of vaccination. Now\n",
    "let's think about the hand-washing campaign. We'll have to answer two\n",
    "questions:\n",
    "\n",
    "1.  How should we incorporate the effect of hand washing in the model?\n",
    "\n",
    "2.  How should we quantify the effect of the money we spend on a\n",
    "    hand-washing campaign?\n",
    "\n",
    "For the sake of simplicity, let's assume that we have data from a\n",
    "similar campaign at another school showing that a well-funded campaign\n",
    "can change student behavior enough to reduce the infection rate by 20%.\n",
    "\n",
    "In terms of the model, hand washing has the effect of reducing `beta`.\n",
    "That's not the only way we could incorporate the effect, but it seems\n",
    "reasonable and it's easy to implement."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "alpine-guest",
   "metadata": {},
   "source": [
    "Now we have to model the relationship between the money we spend and the\n",
    "effectiveness of the campaign. Again, let's suppose we have data from\n",
    "another school that suggests:\n",
    "\n",
    "-   If we spend \\$500 on posters, materials, and staff time, we can\n",
    "    change student behavior in a way that decreases the effective value of `beta` by 10%.\n",
    "\n",
    "-   If we spend \\$1000, the total decrease in `beta` is almost 20%.\n",
    "\n",
    "-   Above \\$1000, additional spending has little additional benefit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "agricultural-trinidad",
   "metadata": {},
   "source": [
    "## Logistic function\n",
    "\n",
    "To model the effect of a hand-washing campaign, I'll use a [generalized logistic function](https://en.wikipedia.org/wiki/Generalised_logistic_function) (GLF), which is a convenient function for modeling curves that have a generally sigmoid shape.  The parameters of the GLF correspond to various features of the curve in a way that makes it easy to find a function that has the shape you want, based on data or background information about the scenario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "blond-armstrong",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import exp\n",
    "\n",
    "def logistic(x, A=0, B=1, C=1, M=0, K=1, Q=1, nu=1):\n",
    "    \"\"\"Computes the generalize logistic function.\n",
    "    \n",
    "    A: controls the lower bound\n",
    "    B: controls the steepness of the transition \n",
    "    C: not all that useful, AFAIK\n",
    "    M: controls the location of the transition\n",
    "    K: controls the upper bound\n",
    "    Q: shift the transition left or right\n",
    "    nu: affects the symmetry of the transition\n",
    "    \n",
    "    returns: float or array\n",
    "    \"\"\"\n",
    "    exponent = -B * (x - M)\n",
    "    denom = C + Q * exp(exponent)\n",
    "    return A + (K-A) / denom ** (1/nu)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "seven-budget",
   "metadata": {},
   "source": [
    "The following array represents the range of possible spending."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "disturbed-amount",
   "metadata": {},
   "outputs": [],
   "source": [
    "spending = linspace(0, 1200, 21)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "supreme-nutrition",
   "metadata": {},
   "source": [
    "`compute_factor` computes the reduction in `beta` for a given level of campaign spending.\n",
    "\n",
    "`M` is chosen so the transition happens around \\$500.\n",
    "\n",
    "`K` is the maximum reduction in `beta`, 20%.\n",
    "\n",
    "`B` is chosen by trial and error to yield a curve that seems feasible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "otherwise-answer",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_factor(spending):\n",
    "    \"\"\"Reduction factor as a function of spending.\n",
    "    \n",
    "    spending: dollars from 0 to 1200\n",
    "    \n",
    "    returns: fractional reduction in beta\n",
    "    \"\"\"\n",
    "    return logistic(spending, M=500, K=0.2, B=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expected-venture",
   "metadata": {},
   "source": [
    "Here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "seventh-strike",
   "metadata": {},
   "outputs": [],
   "source": [
    "percent_reduction = compute_factor(spending) * 100\n",
    "\n",
    "make_series(spending, percent_reduction).plot()\n",
    "\n",
    "decorate(xlabel='Hand-washing campaign spending (USD)',\n",
    "         ylabel='Percent reduction in infection rate',\n",
    "         title='Effect of hand washing on infection rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "legal-michigan",
   "metadata": {},
   "source": [
    "The result is the following function, which\n",
    "takes spending as a parameter and returns `factor`, which is the factor\n",
    "by which `beta` is reduced:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "obvious-congress",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_factor(spending):\n",
    "    return logistic(spending, M=500, K=0.2, B=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sunset-investing",
   "metadata": {},
   "source": [
    "I use `compute_factor` to write `add_hand_washing`, which takes a\n",
    "`System` object and a budget, and modifies `system.beta` to model the\n",
    "effect of hand washing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "polish-multiple",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_hand_washing(system, spending):\n",
    "    factor = compute_factor(spending)\n",
    "    system.beta *= (1 - factor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worthy-fellowship",
   "metadata": {},
   "source": [
    "Now we can sweep a range of values for `spending` and use the simulation\n",
    "to compute the effect:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "statistical-garden",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep_hand_washing(spending_array):\n",
    "    sweep = SweepSeries()\n",
    "    \n",
    "    for spending in spending_array:\n",
    "        system = make_system(beta, gamma)\n",
    "        add_hand_washing(system, spending)\n",
    "        results = run_simulation(system, update_func)\n",
    "        sweep[spending] = calc_total_infected(results, system)\n",
    "        \n",
    "    return sweep"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "clinical-surge",
   "metadata": {},
   "source": [
    "Here's how we run it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "joined-graduation",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import linspace\n",
    "\n",
    "spending_array = linspace(0, 1200, 20)\n",
    "infected_sweep2 = sweep_hand_washing(spending_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "designing-smile",
   "metadata": {},
   "source": [
    "The following figure shows the result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "cheap-decision",
   "metadata": {},
   "outputs": [],
   "source": [
    "infected_sweep2.plot()\n",
    "\n",
    "decorate(xlabel='Hand-washing campaign spending (USD)',\n",
    "         ylabel='Total fraction infected',\n",
    "         title='Effect of hand washing on total infections')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "selective-right",
   "metadata": {},
   "source": [
    "Below \\$200, the campaign has little effect. \n",
    "\n",
    "At \\$800 it has a substantial effect, reducing total infections from more than 45% to about 20%. \n",
    "\n",
    "Above \\$800, the additional benefit is small."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lucky-boulder",
   "metadata": {},
   "source": [
    "## Optimization\n",
    "\n",
    "Let's put it all together. With a fixed budget of \\$1200, we have to\n",
    "decide how many doses of vaccine to buy and how much to spend on the\n",
    "hand-washing campaign.\n",
    "\n",
    "Here are the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "surrounded-copying",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_students = 90\n",
    "budget = 1200\n",
    "price_per_dose = 100\n",
    "max_doses = int(budget / price_per_dose)\n",
    "max_doses"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expired-conditioning",
   "metadata": {},
   "source": [
    "The fraction `budget/price_per_dose` might not be an integer. `int` is a\n",
    "built-in function that converts numbers to integers, rounding down.\n",
    "\n",
    "We'll sweep the range of possible doses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "shaped-utility",
   "metadata": {},
   "outputs": [],
   "source": [
    "dose_array = linrange(max_doses)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "occupational-reply",
   "metadata": {},
   "source": [
    "In this example we call `linrange` with only one argument; it returns a NumPy array with the integers from 0 to `max_doses`, including both.\n",
    "\n",
    "Then we run the simulation for each element of `dose_array`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "recognized-release",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep_doses(dose_array):\n",
    "    sweep = SweepSeries()\n",
    "    \n",
    "    for doses in dose_array:\n",
    "        fraction = doses / num_students\n",
    "        spending = budget - doses * price_per_dose\n",
    "        \n",
    "        system = make_system(beta, gamma)\n",
    "        add_immunization(system, fraction)\n",
    "        add_hand_washing(system, spending)\n",
    "        \n",
    "        results = run_simulation(system, update_func)\n",
    "        sweep[doses] = calc_total_infected(results, system)\n",
    "\n",
    "    return sweep"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cardiac-mitchell",
   "metadata": {},
   "source": [
    "For each number of doses, we compute the fraction of students we can\n",
    "immunize, `fraction` and the remaining budget we can spend on the\n",
    "campaign, `spending`. Then we run the simulation with those quantities\n",
    "and store the number of infections.\n",
    "\n",
    "The following figure shows the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "worth-mounting",
   "metadata": {},
   "outputs": [],
   "source": [
    "infected_sweep3 = sweep_doses(dose_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "adjusted-highlight",
   "metadata": {},
   "outputs": [],
   "source": [
    "infected_sweep3.plot()\n",
    "\n",
    "decorate(xlabel='Doses of vaccine',\n",
    "         ylabel='Total fraction infected',\n",
    "         title='Total infections vs. doses')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dynamic-easter",
   "metadata": {},
   "source": [
    "If we buy no doses of vaccine and spend the entire budget on the campaign, the fraction infected is around 19%. At 4 doses, we have \\$800 left for the campaign, and this is the optimal point that minimizes the number of students who get sick.\n",
    "\n",
    "As we increase the number of doses, we have to cut campaign spending,\n",
    "which turns out to make things worse. But interestingly, when we get\n",
    "above 10 doses, the effect of herd immunity starts to kick in, and the\n",
    "number of sick students goes down again."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "timely-bottom",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose the price of the vaccine drops to $50 per dose.  How does that affect the optimal allocation of the spending?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "swiss-preview",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
